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I. INTRODUCTION 

The slow dynamics of supercooled liquids and glasses 
have been actively discussed for thirty years from the 
view point of both spatial distributions of heterogeneities 
and long time correlations. Recently, the role played by 
heterogeneities in the slow dynamics of glasses has been 
emphasized^ and characterized in terms of 4-point cor- 
relation functions. The long time correlation appear in 
the slow decaying potential part of the shear-stress au- 
tocorrelation function (SACF) and has been called the 
"molasses tail" to differentiate it from the hydrodynamic 
origin of the long time tail in the velocity autocorrelation 
function^ - — and to emphasize its relation to the highly 
viscous glassy stated. The decay of the SACFs have been 
investigated by mode coupling theory (MCT}£ and ki- 
netic theory similar to the velocity autocorrelation func- 
tion that applies to the kinetic part of the shear autocor- 
relation function, but not to the potential part, which is 
central to the molasses tail. Numerical studies eventually 
proved that the long time tail of both the kinetic and po- 
tential parts have a power decay consistent with MCT, 
however, the amplitude of the SACF in dense fluids was 
found to be orders of magnitude greater than predicted. 
This discrepancy was ascribed to the possibility that the 
numerical results were not sufficiently long to resolve the 
long time correlations. However, the cause of this tail is 
likely due to the slow structural relaxation in the dense 
liquid around the peak of the structure factor rather than 
by hydrodynamic phenomena at long wave length. The 
existence of transiently crystal nuclei was also demon- 



a ' Electronic mail: |isobe@nitech ^ac.jp| 
^Electronic mail: [aklerl@llnl.gov] 



stratcd in colloid experiments^ and Langcvin dynamics 
simulation in glassy dense system^. Transient ordering 
mechanism in a quasi-two dimensional liquid near freez- 
ing was also investigated by Sheu and Rice&iii. However, 
the microscopic mechanism of the stress field relaxation 
has not been thoroughly investigated^— . 

Twenty years ago, Ladd and Alder have speculated 
that the long time tail of the shear stress auto-correlation 
function near the solid-fluid transition point in the hard 
sphere system is due to transient crystal nuclei forma- 
tion-^. They found that the potential part of the SACF 
and the angular orientational auto-correlation function 
(OACF) are identical in the long time limit and show 
non- algebraic decay in time. Since the evidence sug- 
gested that the reason for non-algebraic decay is struc- 
tural relaxation rather than hydrodynamic flow, an at- 
tempt was made to understand this slow decaying mech- 
anism by decomposing the OACFs into two-, three-, and 
four-body correlations, however, the four-body correla- 
tions have not been obtained accurately due to computer 
limitation. The prediction of the cooling rate necessary 
to prevent crystallization requires knowledge of the rate 
of growth of a cluster the size of a critical solid nuclei 
and the time they exist in this transient state, and can 
only be obtained from the four body correlation. 

In our previous wor k 16 i , we have investigated the 
slow decay of the pair, Gi, orientational autocorrelation 
function in a two dimensional system consisting of clastic 
hard disks at a single density near the solid-fluid transi- 
tion point, placed in a square box with periodic boundary 
conditions, using a modern fast algorithm based on event- 
driven Molecular Dynamics (MD) simulation^. The 
time evolution of various sized cluster were detected by 
using the bond orientational order parameter— ~—. Near 
the fluid-solid phase transition, we found three regimes 
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in the relaxation of the pair orientational autocorrela- 
tion function, namely the kinetic, molasses (stretched ex- 
ponential), and diffusional power decay (pairs breaking 
apart). We confirmed the non-algebraic decay (stretched 
exponential) at intermediate times presumably due to the 
existence of various sized solid clusters at high densities 
decaying at different rates. 

Then, we focused on the rapidly increasing time with 
increasing density for the decay of the OACFs and were 
able to establish the length of time for which the biggest 
such nuclei exists at each density. The largest cluster near 
the freezing density was found to be only a few sphere 
diameters in size and to persist for typical argon param- 
eters for only about 30 picoseconds. We also compared 
the results to theoretical predictions of the final power 
law decayii. 

To make further quantitative progress, we need to 
investigate the OACF of the quadruplet component, 
C^Ai?, i), as a function of the distance between the 
two colliding pairs AR. From this information it will be 
possible to tell how the cluster size distribution changes 
with time and density, and, subsequently, determine how 
fast one has to increase the density to get a glass in- 
stead of a crystal. Because this is computationally a 
very demanding task, we introduce two methodologies in 
this paper, which are more efficient methods for anal- 
ysis of the quadruplet contribution to the orientational 
auto-correlation function. One is a more efficient coarse 
grained algorithm for calculating pair and quadruplet 
contributions to the OACFs, rather than the previous 
collision based calculation. The other is the extension of 
the usual bond orientational parameter </>g to a higher 
order one involving further neighbor shells. We demon- 
strate that the coarse grained results are in quite good 
agreement with the previous one, but two order of mag- 
nitude faster, allowing for the faster evaluation of the 
4-body autocorrelation functions C^AR, t). 

This paper is organized as follows. In Sec. II, we de- 
scribe how to detect further nearest neighbors and sum- 
marize details of the new algorithm. In Sec. Ill, the re- 
sults of the new improvement are described. Finally, in 
Sec. IV, we summarize the results and discuss the relax- 
ation time of transient clusters and their size distribution. 



II. DETERMINATION OF TRANSIENT CRYSTALS IN 
DENSE LIQUIDS 

In this section, we explain how to categorize neighbor- 
ing shells. Then, generalization of the bond orientational 
order parameter is described in order to calculate the 
autocorrelation of the orientational function. 



A. Detecting Higher Nearest Neighbors 

To consider further neighbors than nearest neighbors 
systematically, we define neighbor shells based on the 
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FIG. 1. Radial Distribution Function at v 
cut-off arrows. 
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1st N.N. 


2nd N.N. 


3rd N.N. 


0.69 


1.579 


2.592 


3.606 


0.65 


1.616 


2.638 


3.660 


0.57 


1.702 


2.747 


3.838 



TABLE I. The distances from the central particle i to further 
neighbor are shown in units of a at various packing fraction. 



minima of the radial distribution functions (RDF) for 
each packing fraction v, which are obtained by an inde- 
pendent calculations via event-driven MD. The radial dis- 
tribution function, g(r), at v = 0.69 in a two-dimensional 
(2D) system composed of N = 4096 hard disks with a di- 
ameter a are shown in Fig. [1] We show 4 red arrows 
for 3 shell radii beyond the central particle, the 1st near- 
est neighbors (N.N.), 2nd N.N., and 3rd N.N., which are 
named by the shell index /, J, K, and L, respectively, 
(e.g., J indicates the particles belong to 1st N.N. against 
the central particle I.) Note that the 2nd N.N. peak has 
a shoulder nearby, which likely indicates that the tran- 
sient crystal becomes significant in dense liquids. The 
result is the shell radii given as the cut-off distance r cut 
for each packing fraction, as summarized in Table |U We 
only consider 3rd N.N. shell in this paper, however, con- 
sidering further neighbors would be straightforward, but 
such large clusters were found to be rare. 

To justify the concept of the above categorization of 
neighbors, we consider the perfect crystal configuration 
for v = 0.69 as a reference, which is shown in Fig. [2] We 
also show 4 red shells (circles) correspond to the cut-off 
radius in Fig. [2] The simulation results on the probabil- 
ity distribution for the particle number of neighbors and 
its averages are shown in Fig.[3]and Table [TTl respectively. 
We hereby justify that neighbors can be detected for 1st 
to 3rd N.N. by setting simple cut-off distances from the 
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FIG. 2. Perfect crystal configuration and neighbor shells at 
v = 0.69. 



V 


1st (J) 2nd (K) 3rd (L) 


perfect crystal 


6 12 18 


0.69 


5.93 11.57 17.23 


0.65 


5.88 11.29 16.71 


0.57 


5.68 10.62 16.32 



TABLE II. The 1st to 3rd N.N. particle numbers for a perfect 
crystal and the actual particle number for various packing 
fractions. 

central particle based on the RDFs. 

B. Decomposition of the Orientational Factors and 
Autocorrelation Functions 

To investigate the temporal properties of clusters in the 
liquid state, the time correlation functions of dynamical 
variables J2i j •^•( r iji*) are considered^ 2 -, where ry in the 
relative distance between the position of two particles 
and Tj. The time correlation function can be written as, 



C(*) = (X;^(ry(0))X;^(rM(<))). (1) 

\ id k,i I 

We have investigated the potential part of the SACF 
(Jxy{t)Jxy(®)) relevant for the molasses tail, where J xy 
is the potential part of the momentum current J% y . For 
a pair- wise potential 4>(rij), J xy is 

where n and m are number density and mass of disks, 
ks and T arc Boltzmann constant and temperature, 




Number of Particles 



FIG. 3. Probability distributions of particle number for each 
neighbor shell at v = 0.69 ( — ) and v — 0.57 (• • • ) are shown. 

( x ijiUij) = { x i ~ x jiUi ~ Uj) are the relative positions 
between particles i and j. In the case of hard disks, this 
becomes, 

£ mb.^MI >-L (3) 

7 

where fry = vy • ry and ^ 7 means the accumulation 
of collisional contributions at the colliding time t 7 . In 
general, C(t) can be decomposed into pair C2(t) (ij — ij 
pair), triplet C^(t) (ij — ik pair), and quadruplet Ci(t) 
(ij — kl pair) contribution a 22 ' 23 , 

c(t) = (x;^«(o))E- A ( r «(*))) w 

\ id M I 

{ift) (kjil) 

= 2J2 <^(ry(0)).A(ry(i))) 

id 

+4 ]T U(rij(0))^(r tt (t))) 

i,j,k 

+ Y, (A(n 3 (0))A(v kl (t))) 

i,j,k,l 

= C 2 (t)+C 3 (t) + C4t). (5) 

Since velocities and positions are no longer correlated 
beyond a few mean collision times, only the orientational 
part of the SACF, namely the orientational autocorrela- 
tion function (OACF), (O xy (t)O xy (0)) , needs to be stud- 
ieSSr—. O xy (t) is defined as 

0^(t) = £^?*(*-M- (6) 

7 
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To avoid the delta function singularity of O xy (t) for hard 
particles, the alternative Einstein-Helfand expression^! 
involving the second derivative, obtained by the numer- 
ical differentiation, is needed for calculating the correla- 
tion function, 



C(t) = (O xy (t)O xy (0)) 



(7) 



needs only configurational data of particles at a certain 
discrete time t = A x At instead of collisions, where A is 
an integer number and At is taken at an interval e.g. the 
mean free time, At = to. The essential idea of our im- 
provement is the coarse graining (CG) of collision events 
in neighbor shells at discrete small times, in which we 
recognize the tagged particles as the candidates of col- 
lisions within At. Therefore, we can calculate orienta- 
tional factors by summing those candidates with the fol- 
lowing slight modification, 



where 



G(t) 



(8) 



and where O(i) is the unit step function. Note that there 
are three independent oricntational factors {O xy , yzi 
and Ozx) in 3D. The pair and quadruplet contributions 
of OACF arc defined as, 



C 2 (t) 



2dt 2 



1 d 



N 



\ 1,3 
I N 



Ci{t) ~2d£\ £ G ij (t)G kl (t) 



(9) 



(10) 



X y(t) = £ 



xu{t)yu(t) 
^ (<r + 5r) 2 

£cos0/j sin Qjji 
i.j 



(12) 



where Sr = \/xu(t) 2 + yu(t) 2 — a (<C a) is the small 
gap distance between particle / and J and 6u is the an- 
gle of vector r/j against a reference axis (e.g. x-axis). 
Such a modification has a great advantage in efficiency 
since we do not wait for the actual particle collisions as 
in the event-driven scheme. We restrict the particle pairs 
considered as the reference pairs at the start of each simu- 
lation as nearest neighbors. This is expected to improve 
sampling drastically by quick sorting over choosing ar- 
bitrary reference particle pairs as we have done in the 
previous methodic. 



since (G y '(0)) = (G fei (0)) = (G lJ (t) = 0%(t),G kl (t) = 
O x l y (t)). To ease calculating C 2 and C4, we introduce a 
"collision pair index" 



lk = ( 7i - l)N - 7<( 7i - l)/2 + 7j - 7i 



(11) 



where ji and jj are particle indexes of colliding pairs, 
which identifies a given pair quickly, thus avoiding hav- 
ing to check whether the same collision pair has collided 
before. For example, in case of JV = 4, the total num- 
ber of collision pairs are N(N — l)/2 = 6, which can be 
listed as ( 7< , 7i ) = (1, 2), (1,3), (1, 4), (2, 3), (2, 4), (3, 4), 
where 7 j < jj . By using the "collision pair index" , we 
obtain 7 fc = 1,2,3,4,5,6 for each collision pair, respec- 
tively. Therefore, it is convenient to deal with the colli- 
sion pair as the sequential number to sort and insert into 
the array of correlation pairs of G l: >(t). This speeds up 
the calculation considerably. All properties in this paper 
are normalized and indicated by C*(t*), where t* is the 
reduced time t* = t/to(to is the mean free time). 



C. Higher Order Orientational Factor Based on Course 
Graining Method 

Here, we propose an alternative methodology for cal- 
culating the OACF efficiently, which differs from the pre- 
vious method based on collision events. The new method 



D. Generalized Order Parameter Based on Crystal 
Structure 



1. Bond Orientational Order Parameter 

The usual bond orientational order parameter 4>6 for a 
hard disk i is defined by^, 



1 Nt 

— Eexp(6i6 



(13) 



where Ni is the number of the nearest neighbors around 
the tagged particles i, and Q\ is the angle between the 
position vector from the disk j to i and an arbitrary fixed 
reference axis (e.g., x-axis). 

Previously, two disks are defined as nearest neighbors if 
the separation is within a distance, say 1.4 to 1.7a. This 
choice is reasonable from the view point of our definition 

by the RDFs. The absolute value $| = 



; 0g takes on 

values between and 1 , and measures the degree of crys- 
tallization in terms of considering only nearest neighbors, 
where (j>\ is the complex conjugate of 4> l 6 . Figure 0] shows 
typical snapshots of the spatial distribution of $| at the 
packing fraction v = 0.65 (left) and 0.69 (right), respec- 
tively. The gradation in shading of the particles indicates 
the value of $g; the darker, the closer to unity. We clearly 
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0.5 



FIG. 4. The spatial distribution of $6 for 4096 particles system at a given time for two packing fractions, v = 0.57(left) and 
0.69(right). The darker the region, the closer $6 is to unity. 



central particle (I) (0,0) 

1st N.N. (J) (1,0), (0,1), (-1,1), (-1,0),(0,-1), (1,-1) 
2nd N.N. (K) (1,1), (-1,2), (-2,1), (-1,-1), (1,-2), (2,-1) 

(2,0), (0,2), (-2,2), (-2,0), (0,-2), (2,-2) 

3nd N.N. (L) (3,0), (2,1), (1,2), (0,3), (-1,3), (-2,3) 

(-3,3), (-3,2), (-3,1), (-3,0), (-2,-1), (-1,-2) 

(0,-3), (1,-3), (2,-3), (3,-3), (3,-2), (3,-1) 



TABLE III. (C, r?) pairs for neighbors. 

observe the dramatic growth of several solid nuclei as the 
density nears solidification. However, those solid nuclei 
will disappear after a certain transient timeii. To in- 
vestigate those transient clusters more quantitatively, we 
extended the usual bond orientational order parameter 
toward further neighbors, that is, the 2nd to 3rd N.N. 
and investigate the relaxation time via autocorrelation 
functions. 



2. Angle for Neighbors based on Triangle Crystal Lattice 

Here, the angle between bond vectors for the general- 
ized order parameter is considered. It is easy to deal with 
the neighbors when we introduce the integer index (£, rf) 
summarized in Table Mil If the particle is located on the 
crystal lattice with the side distance a, the position vec- 
tor of particle c can be described as 

c = Ca + rib, (14) 

where (C,v) are integers and a = (l,0)a, b = (|, y ^)a 
are unit vectors of the crystal lattice. 



We next consider the angle of bond vectors. If the 
particle positions are in a perfect crystal, that is r = c, 
the bond vector between central particle I and 1st N.N. 
J can be defined as, 

cji = cj-a. (15) 

We can define 6 kind of bond vectors (cj/, cki, cli ckj, 
Clj, and c^k) within 3rd N.N. Angles 9 between bond 
vectors for cjj and cj/j (J' 7^ J) are easily calculated as 
8 = n' x 7r/3, where n' is an integer, which reduces to the 
usual calculation for <j)Q. In general, the angle between 
bond vectors can be calculated by, 

We can define 432 pairs of bond vectors within 3rd 
N.N., of which 6 are for cjj, 12 for cki, 18 for cli, 72 
for ckj, 108 for clj, and 216 for c^k- We don't consider 
ckj, clj in this paper. Based on the perfect crystal, we 
calculate the angle probability distribution between the 
bond vector pairs (i) (cji, Cji) (c K i, c K i) (c LI , c LI ), 
(ii) (cjj, ckl)- Thus, the concept of bond orientational 
order parameter for further neighbors can be extended. 

3. Generalized Order Parameter Based on Crystal 
Structure 

The usual 4>e order parameter is obtained from the 
angle between the position vector from disk J to I and 
"an arbitrary fixed reference axis" . To consider the an- 
gle between 1st and 2nd neighbors pairs (and more), the 
bond angle can be redefined by relative vectors between 
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1-3 and I-K (etc..), instead of "an arbitrary fixed refer- 
ence axis" . The complex generalized order parameter for 
a tagged particle I with the actual position vectors rjj 
can be generalized by 



ti= NjiNj-l) Ex.(rjJ,r^j), (17) 
X s (rj7,rj-/) =exp(isO(rji,Tjii)), (18) 
e(Tj Z ,rj fI ) = 9 J / = (9{ 6j) = cos- 1 ^ gLg|^ >) 

If s = 6 and tcjii = (I, 0) (i.e., unit vector of x-axis) are 
fixed, we can reduce the above expression to the usual (fie 
order parameter described in eq. (|T3|) . 



III. RESULTS 

The systems considered are hard disks placed in a 
L x x L y (= A) square box with periodic boundary condi- 
tions. Initially, the simulation systems for each packing 
fraction v (— Nir(a /2) 2 /A) are prepared in an equilib- 
rium state by a sufficiently long preliminary run. For the 
close packed area A , v = n/(2\/3(A/Ao)). The system 
evolves through collisions, using an algorithm based on 
event-driven MD simulation-^. Most of the calculations 
are done with a relatively small particle numbers of par- 
ticles, N = 4, 096, since long time runs for calculating ac- 
curate tails of autocorrelation functions are needed. We 
confirmed that periodic boundary effects on the OACF 
don't appear for this system sizei£. The density is set 
at relatively dense values near the solid-fluid transition 
point v c near 0.70, namely at v = 0.69, 0.65, and 0.57, 
primarily to compare with the previous result o 16 ' . 



A. Total and Pair Orientational Autocorrelation Functions 

The algorithm for C2{t) is described in the following 
steps, 

f . Prepare the list of vectors for particle pairs within 
the 3rd N.N. shells at t = 0, which are likely col- 
lision candidates in a short time. This is called "a 
reference particle pair list". This list is similar to 
the neighbor list known in the standard MD tech- 
nique to increase efficiency. Note that if we only 
register particles within I N.N. shell, they can es- 
cape from that shell after a short time and there- 
fore, the 3rd N.N. choice is a better one. 

2. Sort the above list vectors by sequential number 
according to the collision pair index 7^. (eq. (fTTj)). 
where 7^ and jj are the index of the colliding can- 
didate particles. 




t* (-Mb) 

FIG. 5. Comparison on the OACF (C to tai) at two densities 
between previous collision-based method and new method. 
The vertical axis is normalized and the horizontal axis is the 
scaled time by the mean free time, to- 



3. During the event-driven MD simulation at each dis- 
crete time t = AAi (A = 1, 2, 3, • • • ), we list particle 
pairs within 3rd N.N. shell into such arrays and sort 
them sequentially by using the pair list index. 

4. Next, compare the pair particle index 7^ obtained 
from configurations at time t with the list of refer- 
ence pair index prepared at time t = 0. 

5. If the pair particle index is found to be the same as 
the reference pair index for the 1st N.N., we insert 
it in the CG modified orientational factor O l J y (t) 
(eq. (fT2"|) b If the particle index is not the same as 
in the list of reference pairs, we discard it. 

6. After the simulation is performed for a long time 
we obtain averages for 7^ and hence C2{t) via the 
Einstcin-Hclfand formula with time resolution At. 

In Figs [5] and [6l comparison of the C total and C2 at two 
densities between the previous collision-based method 
and the new method are shown. The relaxation time 
and long time behavior of both Ctotai and C2 obtained 
by the new method arc in fairly good agreement with 
that of the previous method at both densities^. How- 
ever, the efficiency proved to be drastically improved. 
Although the efficiency depends on the particle number 
N which determines the total array for particle pairs, the 
efficiency was more than 16 and 77 times faster than pre- 
viously for the total and pair autocorrelation function at 
[N, v) = (4096,0.69), respectively. 
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t* (-Mb) 

FIG. 6. Comparison of C2 at two densities between previous 
collision-based method and new method. The vertical axis is 
normalized and the horizontal axis is the scaled time by the 
mean free time, to- 



B. Quadruplet Orientational Autocorrelation Functions 

We then showed that the new method can calculate 
higher order correlations such as the distance dependence 
of the autocorrelation for the quadruplet contribution, 
which cannot be resolved by the collision based calcula- 
tion^. 

The C4 algorithm is basically constructed as was C2, 
however, it is more complex since we need to use more 
sorting procedures and searching method for detecting 
valid pairs for i — j and k — I within 3rd N.N. The main 
computational task is to search the particle pair index 
in the reference particle pair list for both i — j and k — 
I. To detect the valid particle pairs for i — j and k — I 
within 3rd N.N. separated by a calculable distance at 
each discrete time is relatively easy, since those candidate 
pairs have already been registered at t — 0. Therefore, it 
leads to a significant improvement for calculating the C4 
autocorrelation functions. Note that since not only i — j 
pair but also k — l must be the collisional candidates, k — l 
pairs are searched under the condition that k — I are also 
1 N.N. 

Figure [7] shows the probability distribution function of 
AR in units of a, where AR is the separation between 
the center of mass of the particle pairs i—j and k — l. The 
contribution beyond the 4th N.N. are not calculated. For 
a perfect crystal, the number of distinct quadruplet pairs 
(i — j,k — I) for tagged particle i are ~ 118. 5iV in the 
N particle system. The actual sampling of quadruplet 
pairs are about 91 pairs for each tagged particle i at v = 
0.69. Figure [5] shows the time dependence of C^AR^t) 
normalized by Ci(AR, 0) for two packing fraction v = 
0.69 and 0.65 for the AR/a corresponding to 1st, 2nd, 



T 




' ' ' 2 ' ' ' 4 
AR ([a]) 



FIG. 7. The probability distribution function of AR in the 
unit of a for v = 0.65 and 0.69. 



and 3rd peaks in Fig. [7J respectively. 

In Ci{AR, t) beyond AR/a ~ 3.4, the contribution of 
outer particles of 4rd N.N. shell gradually becomes dom- 
inant. Computational costs are just a few times larger 
than for the C2. In Fig. [5J we found that C4 for 1st 
peak changes from positive to negative values around 
t* ~ 246(^ = 0.69) and 53(^ = 0.65), respectively. On 
the contrary, C4 for 2st peaks changes negative to posi- 
tive values around t* ~ 197(V = 0.69) and 39(^ = 0.65), 
respectively. This is because the geometry of configura- 
tions between i — j and k — I pairs. Those results indi- 
cate that C4 in v = 0.65 decays much faster than that 
of v = 0.69. In our new method, Cn{AR,t) can resolve 
how the cluster size distribution changes in time for each 
density quantitatively. 



C. Generalized Order Parameter 

Tables IIVI shows the results for the generalized order 
parameters (GOP), = V^> l s , for s = 6, 12, 18 (i.e., 
1st -3rd N.N.) as described in Sec. II. D for each packing 
fraction v = 0.57, 0.65, and 0.69. In Table El the GOPs 
are calculated based on the fixed reference axis (i.e. x- 
axis, Tj'i = (1,0)). We note that all GOPs decrease for 
higher neighbors for each packing fraction and increase 
for the higher packing fraction. Figure [9] shows the spa- 
tial distribution of GOP $ l s (s = 12 (left) and s = 18 
(right)) for 4096 particles system at a given time at the 
packing fractions v = 0.69. Comparing with Fig. @J the 
darker the region in Fig. [4] gradually decreases for 2nd 
and 3rd N.N. This can be useful to determine the crys- 
tal size quantitatively and to further characterizing the 
transient clusters. 
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t* (-Mb) 

FIG. 8. The distance dependence normalized d(AR, i) for 
each packing fraction v = 0.69 and 0.65 are shown. 
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0.57 


0.65 


0.69 


<<2>6> 


0.487 0.556 


0.645 


<-3>! 2 ) 


0.369 


0.416 


0.496 


(<J>18> 


0.300 


0.341 


0.402 



TABLE IV. The generalized order parameter, at various 
packing fractions. 



D. Autocorrelation Function Based on Generalized Order 
Parameter 

The autocorrelation functions based on the generalized 
order parameter for each particle i is decomposed into 
real and imaginary parts, 

#(*) = Re(#(*))+am(#(*)), (20) 
leading to the following autocorrelations, 



C GOP (t) - (CW s (0)) (21) 
= (Re(^ s (t))Re(^(0))+Im(^(i))Im(^(0^) 

The decay of this autocorrelation function CqopM 
for s = 6, 12, 18 at v = 0.69 and 0.65 normalized by 
C GOP (t) = C GO p(t)/(® l s ) 2 is given in Table El The de- 
cay of the shear stress is directly related to the decay of 
C (see, Fig. |5]and Table W\ . and its life time to the disso- 
lution of the cluster, i.e., when its core, C2, melts 16 d 7 . A 
movie of a transient crystal nuclei formation shows that 
the growing process involves the increasing order of the 
neighboring particles to as large as 3rd neighbors followed 
by the dissolving process till the nearest neighbors are no 



V 


T{C t otal) 


r(C 2 ) 


n(C 4 ) 


T 2 (C 4 ) 


r GOP 


•2 

r GOP 


T GOP 


0.69 


12 


133 


78 


66 


40 


20 


15 


0.65 


6 


41 


21 


21 


9 


7 


6 



TABLE V. Relaxation time r for C to tai, C2, and Tk of 
C(AR, t) for 1st (k = 1) and 2nd (k = 2) peaks, and tqop 
for n— th N.N. at v — 0.69 and 0.65. All time are in the unit 
of to (mean free time). 

longer ordered. The overall average of all nucleation pro- 
cesses has a time scale given in Table fVl at v = 0.69 of 
40, 20, 15 for the 1st, 2nd, and 3rd N.N., respectively. 

IV. SUMMARY AND DISCUSSION 

In this paper, a method for analysing higher order pa- 
rameter of the liquid state is developed, especially to in- 
vestigate transient crystals in dense liquid systems. In- 
stead of calculating orientational factors and their de- 
composition based on collision event as previousl y 16 ' 17 , 
we developed a more efficient methodology by introduc- 
ing neighbor shells which coarse grain collision events 
within these neighbor shells as candidates of further col- 
lisions during short times. We confirmed that the same 
results and relaxation time are obtained as the previ- 
ously. We also demonstrate that this improvement per- 
mits calculating higher order parameters of orientational 
contributions such as the distance dependence of the au- 
tocorrelation function of the quadruplet contributions 
Ci(AR, t), providing information on the size and num- 
ber of the transient crystals and their life time. The 
size and number distribution can be investigated not only 
for 1st N.N. but also for further neighbors, as shown in 
Fig. 1101 If we somewhat arbitrary recognize particles 
with > 0.9 as the center of crystal-like structures as 
before^, it is possible to estimate the crystalline fraction 
for higher neighbors, as shown in Fig. 1111 

At v = 0.69, we found that ~ 15% of the particles in 
the system have crystal-like nearest neighbors, ~ 1.8% 
2nd N.N. and ~ 0.13% 3rd N.N. This means that there 
are ~ 5 (i.e., = 4096 x 0.0013) crystal clusters of the 
size of 3rd N.N. shell or the size of ~ 36 particles in the 
system. The typical snapshot of the spatial distribution 
confirmed the above estimation of size and cluster num- 
ber. (See, the right of Fig. [TT]). Rarely have the central 
particles and some of the next nearest neighbor order be- 
yond > 0.9 at the same time because such 3rd order 
N.N. are rare events. 

The lifetime of these clusters given in the unit of to 
(mean free time) are determined from the autocorrela- 
tions for C to tai{t), C2(t), C4(AR, t) and Cgop(G that are 
of the stretched exponential form. The relaxation time 
is defined as the time when the auto-correlation function 
decays to 1/e of its initial value. In Table [V] the relax- 
ation time r for Ctotai, C2, and of C(AR,t) for k-th 
peaks, and r GO p for n-th N.N. at v = 0.69 and 0.65 are 
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FIG. 9. The spatial distribution of generalized order parameter &l (for s = 12 (left) and s = 18 (right)) in 4096 particles 
system at a given time at v = 0.69. The darker the region, the closer is to unity. For s = 6, see Fig. fright). Note that 
GOPs of higher order N.N. for a given central particles may occasionally have a higher correlation than that of a lower order 
N.N. because the core particles may not be as highly ordered. 




second peaks although the amplitudes arc different, while 
Ci(AR : t) for the 1st peak at v = 0.69 has a somewhat 
slower decay than that for the 2nd peak. The rapid in- 
crease in relaxation time of clusters with density as well 
their number and size is closely related to the rapid in- 
crease in viscosity near freezing. If we cool the system 
or compress much faster than the relaxation time esti- 
mated by the methods presented here, we can determine 
the condition under which glass might form. 



GOP H> 8 '| 



FIG. 10. The probability distribution functions of &l at 
v = 0.69 for further neighbors. 



summarized. As expected, the relaxation time increase 
when the packing fraction increases and decreases for the 
higher neighbor shells. In comparing the relaxation time 
for C2, C4 for the 1st peak and GOP for 1st N.N., ac- 
count must be taken of the number of particles involved, 
namely 2, 4 and 6 particles respectively. The larger the 
numbers, the quicker their order is destroyed. C^AR, t) 
for the 1st peak loses orientational order faster than C2 
and Tqqp even faster. Interestingly the relaxation curves 
for C4 at v — 0.65 are almost the same for the first and 
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FIG. 11. (left) The percentage of GOP with $ l s > 0.9 in terms oin-th N.N. for v = 0.65,0.69 are shown, (right) The 
distribution of GOP for the central particle of a cluster is plotted using threshold V&* > 0.9 for 3rd N.N. in the same configuration 
as Fig. 9. 
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